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1  Introduction 


The  study  of  path  effects  of  complex  structure  and  heterogeneities  on  the  excitation 
and  propagation  of  regional  phases  in  different  areas  remains  critical  for  both 
discrimination  and  yield  estimation  procedures  for  monitoring  the  CTBT.  The 
problem  will  be  most  severe  in  the  case  of  Non-Proliferation  monitoring,  in  which 
the  potential  nuclear  tests  may  occur  in  very  different  geological  and  geophysical 
environments.  Today,  regional  waves  are  one  of  the  most  important  information 
sources  for  monitoring  purpose.  Due  to  the  complexity  involved  in  regional  phase 
propagation,  synthetic  simulations  will  play  an  important  role  in  areas  where  there 
is  a  lack  of  sufficient  observations.  To  meet  these  requirements,  the  ultimate  goal 
is  to  develop  a  computationally  viable  technique  for  calculating  high-frequency 
(1  -  25  Hz)  synthetic  seismograms  in  regional  distance  (  >  1000  km)  for  three- 
dimensional,  heterogeneous  (on  large  and  small  scales)  crustal  structures  including 
rough  surface  and  interfaces. 

In  the  past,  boundary  integral  equation  (BIE)  or  boundary  element  (BE)  meth¬ 
ods  have  been  extensively  used  to  study  the  effects  of  topography  or  sedimentary 
basin  structures  on  ground  motions  at  the  surface.  These  have  also  been  used  to 
study  the  Lg  blockage  problem  with  limited  success.  Blockage  is  assumed  to  be 
caused  by  coastlines,  mountains  and  sudden  change  of  crustal  thickness.  However, 
numerical  simulations  of  blockage  by  large-scale  crustal  structures  have  not  suc¬ 
ceeded  in  matching  the  observations  (Campillo  et  al.,  1993;  Gibson  and  Campillo, 
1994).  Most  simulations  are  either  for  surface  topography  or  for  irregular  structure 
beneath  a  flat  surface  (sedimentary  layer)  due  to  the  restriction  of  computational 
complexity.  However,  the  combination  of  both  surface-topography  and  sedimen¬ 
tary  structure  may  have  more  dramatic  influence.  An  irregular  surface  and  low- 
velocity  layer  can  both  trap  part  of  the  Lg  energy  into  the  surface  layer  and  scatter 
the  Lg  wave  out  of  the  crustal  waveguide.  Existing  methods  are  also  not  capable 
of  simulating  the  combined  effects  of  both  large-scale  structure  and  the  associ¬ 
ated  small-scale  heterogeneities.  Irregular  topography  and  near-surface  structure 
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are  the  manifestation  of  past  and/ or  present  tectonic  processes  which  often  pro¬ 
duce  crustal  heterogeneities  at  different  scales.  The  effects  from  the  small-scale 
(wavelength-scale)  heterogeneities  must  be  taken  into  consideration  in  modeling 
blockage  and  other  Lg  propagation,  scattering  and  attenuation  phenomena. 

Our  aim  is  to  develop  a  new  hybrid  numerical  method  by  combining  the  general¬ 
ized  screen  method  with  the  boundary  integral  equation  method.  The  generalized 
screen  method  can  handle  wave  propagation  in  a  heterogeneous  waveguide  with 
modest  topography.  The  method  is  based  on  one-way  wave  equation  theory  (Wu, 
1994;  1996;  Wu  and  Xie,  1994;  Wu  and  Huang,  1995).  In  the  crustal  waveguide 
environment,  major  wave  energy  is  carried  by  forward  propagating  waves,  includ¬ 
ing  forward  scattered  waves,  and  therefore  the  neglect  of  backscattered  waves  in 
the  modeling  will  not  change  the  main  features  of  regional  phases  in  most  cases. 
By  neglecting  backscattering  in  the  theory,  the  method  becomes  a  forward  march¬ 
ing  algorithm.  The  present  value  of  the  wavefield  in  a  vertical  cross-section  in 
the  waveguide  can  uniquely  determine  the  wavefield  in  the  next  vertical  cross- 
section.  In  this  way,  we  can  obtain  the  wavefield  in  successive  vertical  cross-sections 
throughout  the  entire  waveguide.  This  method  provide  us  with  enormous  savings 
in  computing  time  and  storage,  and  makes  it  a  very  efficient  method  which  can 
propagate  high  frequency  regional  signals  to  very  long  distances. 

Modest  surface  topography  can  be  modeled  by  coordinate  transformation  in  the 
generalized  screen  method.  However,  the  algorithm  for  handling  the  topography 
is  still  in  the  process  of  development.  On  the  other  hand,  the  boundary  integral 
equation  and  boundary  element  methods  have  the  flexibility  needed  to  incorporate 
complex  topographic  features  into  the  model.  However,  since  matrix  operations 
are  involved,  the  boundary  integral  equation  method  is  not  that  efficient.  When 
the  ratio  of  model  dimension  to  wavelength  is  too  large,  the  computation  time 
and  memory  requirement  become  formidable.  This  problem  can  be  circumvented 
through  a  hybrid  method.  The  hybrid  method  will  combine  the  advantages  of 
the  above  mentioned  two  methods  and  avoid  their  disadvantages.  The  Lg  phases 
generated  by  the  source  are  propagated  to  a  certain  distance  with  the  generalized 
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screen  method.  Then,  the  output  will  be  used  as  the  input  to  boundary  integral 
equation  (BIE)  or  boundary  element  (BE)  methods,  and  the  later  is  used  to  cal¬ 
culate  the  interaction  between  the  Lg  wave  and  the  complex  waveguide  structure 
with  rough  topographic  features.  This  approach  provides  the  ability  to  investigate 
the  interaction  between  the  Lg  wave  and  crustal  waveguides  having  complicated 
structures  including  severe  topography  for  long  distance  propagation. 

In  this  project,  the  screen  method  has  been  successfully  developed  for  a  crustal 
waveguide  for  2-D  SH-wave  propagation.  The  boundary  integral  equation  and 
boundary  element  methods  and  the  corresponding  connection  schemes  to  the  screen 
propagator  have  been  developed  and  tested.  Numerical  examples  demonstrated 
the  feasibility  of  this  hybrid  approach.  Preliminary  numerical  simulations  on  Lg 
propagation  through  complex  waveguides  with  rough  surface  or  rough  bottom  of 
sedimentary  layers  show  interesting  results. 

2  Generalized  Screen  Method 

For  an  isotropic  2D  elastic  medium,  the  SH  and  the  P-SV  waves  are  decoupled. 
Here,  we  treat  only  the  SH  problem  to  demonstrate  the  applicability  of  the  screen 
propagators  to  a  crustal  waveguide.  Under  such  a  circumstance,  the  equation  of 
motion  becomes 


-o,V(rWr)  =  (1) 

where  to  is  the  frequency,  r  =  ( x ,  z)  is  a  2D  position  vector,  u  is  transverse  displace¬ 
ment,  p  is  the  density  of  the  medium,  and  p  is  the  shear  rigidity.  We  decompose 
the  parameters  of  the  elastic  medium  and  the  total  wave  field  into 

P  =  Po  +  tip 

P  =  Po  + 

u  =  u°  +  U  (2) 
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where  p0  and  p0  are  parameters  of  the  background  medium,  8p  and  Sp  are  corre¬ 
sponding  perturbations,  u°  is  the  primary  field  and  U  is  the  scattered  field.  Then 
the  SH  wave  equation  can  be  rewritten  as 


Mo V2lf  +  u2p0U  =  ~[u28pu  +  V  •  SpVu] , 

(3) 

or 

(V2  +  k2)U(r)  =  -k2F(r)u(r) , 

(4) 

where  k  =  ujv  is  the  wavenumber  in  the  background  medium  and  v 
ground  S  wave  velocity  defined  by 

is  the  back- 

C5 

II 

0 

0 

(5) 

In  the  right-hand  side  of  (4),  F(r)  is  a  perturbation  operator 

E(r)  =  £,(r)  +  ^V-eMV, 

(6) 

with 

e?(r)  =  , 

Po 

(7) 

,  .  Sp( r) 

£M(r)  =  -Z±-L  . 

Mo 

(8) 

Eq.  (4)  is  a  scalar  Helmholtz  equation.  With  a  half-space  scalar  Green’s 
function  gh,  the  scattered  field  U  can  be  written  as 


^(ri)  =  k*  jv  d2vgh{rl-r)F(r)u{r) ,  (9) 

where  the  2D  integration  is  over  the  volume  V  including  all  the  heterogeneities 
m  the  modeling  space.  Under  the  forward-scattering  approximation,  the  total 
field  and  Green’s  function  under  the  integration  in  the  above  equation  can  be 
replaced  by  their  forward-scattering  approximated  counterparts,  and  the  field  can 
be  calculated  by  a  one-way  marching  algorithm  along  the  x-direction  using  a  dual 
domain  technique. 


4 


2.1  Wide-Angle  Screen  Approximation 

The  half-space  model  can  be  sliced  into  thin-slabs  perpendicular  to  the  propagation 
direction.  The  weak  scattering  condition  holds  for  each  thin-slab.  For  each  forward 
step,  the  forward-scattered  field  by  the  thin-slab  is  calculated  and  added  to  the 
primary  field  so  that  the  updated  field  becomes  the  incident  field  for  the  next 
thin-slab.  The  formulas  of  the  dual-domain  implementation  are  summarized  as 
follows: 


U(xuKz)  =  U„(xi,Kz)  +  U^xuKz) 


where 


Up(xi,  I\z)  =  ik  f 1  dxe^x^C[-£p(z)u0{z)] 

Jx'  7 

Kz)  =  ik  r  dxe^-^  {c[e/i(z)a*u0(^)]  -  iS[— e^z)Bzu0{z)] 
Jx'  {  7 

where  C[f(z )]  and  <S[/(.j)]  are  the  cosine  and  sine  transforms,  defined  by 

roo 

C[f(z)]  =  /  dz2  cos(K2z)f(z) 

Jo 

roo 

S[/(z)]  =  /  dz2sin(Kzz)f(z) 

Jo 

In  Eq.  (11)  and  (12),  u0,  Bxu0  and  Bzuo  can  be  calculated  by 

1  r  C1© 

uo(x,z)  =  —  /  dK'eiK:ze{^x-xr>u0(x',K') 

Z7T  J — oo 

=  C-^e^-^Uoix'jC)} 


and 


Bxu0(x,z)  =  C-i[e’V(*-*')Xw  0{x'jC)] 
Bzu0(x,z)  =  iS-1[e*(*-x')j±uo(x,,K')] 


(10) 


(11) 

(12) 


(13) 


(14) 


(15) 
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Eq-  (11),  (12),  (14)  and  (15)  are  the  dual-domain  expressions  of  the  wide-angle 
screen  propagator  for  half-space  SH  problems. 

The  procedure  can  be  summarized  as  follows. 

1.  Cosine  transform  the  incident  fields  at  the  entrance  of  each  thin-slab  into 
the  wavenumber  domain. 

2.  Free  propagate  in  the  wavenumber  domain  and  calculate  the  primary  field 
and  its  gradient  within  the  slab. 

3.  At  each  horizontal  position  within  the  slab,  inverse  cosine/sine  transform 
the  primary  field  and  its  gradients  into  the  space  domain  and  then  interact  with 
the  medium  perturbations  ep  and  ep. 

4.  Cosine/sine  transform  the  distorted  fields  into  the  wavenumber  domain  and 
perform  the  divergence  operations  to  get  the  scattered  fields. 

5.  Calculate  the  primary  field  at  the  slab  exit  and  add  to  the  scattered  field  to 
form  the  total  field  as  the  incident  field  at  the  entrance  of  the  next  thin-slab. 

6.  Continue  the  procedure  iteratively. 

2.2  Small-Angle  Screen  Approximation  and  the  Phase-Screen 
Propagator 

When  the  energy  of  crustal  guided  waves  is  carried  mainly  by  small-angle  waves 
(with  respect  to  the  horizontal  direction),  small  angle  approximations  can  be  in¬ 
voked  to  simplify  the  theory  and  calculations.  Under  the  phase-screen  approx¬ 
imation,  the  heterogeneous  half-space  is  represented  by  a  series  of  half-screens 
embedded  in  the  homogeneous  background  half-space.  Waves  propagate  between 
screens  in  the  wavenumber  domain  and  interact  with  phase-screens  in  the  space 
domain.  The  interaction  is  only  a  phase-delay  operator  (multiplication  in  space 
domain).  The  formula  for  dual-domain  implementation  is 

u(xuKz)  =  u0(xuK-)  +  U(xuKz) 

=  e h{xi~x,)  I"  dz2cos(Kzz)[l  +  zk2Ss(Z)]u0(x',  z) 

J  0 
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0il{xi-x‘ 


,2  %kS,(z) 


Uo{x  , 


(16) 


where  exp[2iA:Ss(z)]  is  the  phase  delay  operator.  The  procedure  can  be  summarized 
as  follows. 

1.  Cosine  transform  the  incident  field  at  the  starting  plane  into  the  wavenumber 
domain  and  free  propagate  to  the  screen. 

2.  Inverse  cosine  transform  the  incident  field  into  the  space  domain  and  interact 
with  the  shear  slowness  screen  (phase-screen)  to  get  the  transmitted  field. 

3.  Cosine  transform  the  transmitted  field  into  the  wavenumber  domain  and 
free  propagate  to  the  next  screen. 

4.  Repeat  the  propagation  and  interaction  screen-by-screen  to  the  boundary  of 
the  model  space. 


2.3  Treatment  of  the  Moho  Discontinuity 

The  Moho  discontinuity  can  be  treated  in  two  ways.  One  is  to  put  the  impedance 
boundary  conditions  in  the  formulation,  the  other  is  to  treat  the  parameter  changes 
as  perturbations  which  are  therefore  incorporated  into  the  screen  interaction.  The 
former  has  the  advantage  of  computational  efficiency.  The  latter  has  the  flexibility 
of  handling  irregular  interfaces.  Here,  we  adopt  the  latter  approach  and  check 
the  validity  of  perturbation  approach  for  the  Moho  discontinuity  by  a  reflectivity 
method  and  a  finite  difference  algorithm. 

2.4  Numerical  Tests  for  the  Screen  Method 

In  this  section  we  give  examples  of  using  the  half-space  phase-screen  algorithm 
for  regional  wave  propagation.  First,  we  show  the  accuracy  of  the  method  by- 
comparing  the  synthetic  seismograms  generated  by  the  screen  method  with  those 
calculated  by  the  reflectivity  method.  Shown  in  Figure  1  are  synthetic  seismograms 
for  a  simple  one  layer  crust  model  (a  ID  model)  with  a  flat  Moho  discontinuity 
at  the  depth  of  32  km.  The  source  function  is  a  Ricker  wavelet  with  a  dominant 
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Figure  1:  Comparison  of  synthetic  seismograms  along  the  surface  calculated  by 
the  screen  method  and  reflectivity  method  in  wave-number  domain.  The  model  is 
a  simple  one  layer  crust  model  with  a  flat  Moho  discontinuity  at  the  depth  of  32 
km.  The  source  function  is  a  Ricker  wavelet  with  a  dominant  frequency  of  1.0  Hz. 
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SP(thick  solid  line)  &  FD(thin  line)  for  inhomogeneous  model 


Figure  2:  Comparison  of  synthetic  seismograms  along  a  vertical  profile  at  a  dis¬ 
tance  of  250  km  calculated  by  the  screen  method  and  a  finite-difference  method. 
Shown  in  the  upper  panel  is  the  crustal  model  and  the  lower  panel  shows  synthetic 
seismograms.  The  thin  lines  are  from  the  finite  difference  method  and  the  thick 
lines  are  from  the  generalized  screen  method. 
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frequency  of  1.0  Hz.  Thin  lines  are  from  the  reflectivity  method  and  thick  lines 
are  from  the  generalized  screen  method.  The  source  is  located  at  a  depth  of  2 
km,  and  its  time  function  is  a  Ricker  wavelet  with  a  dominant  frequency  of  1.0  Hz. 
Except  for  near  vertical  reflections,  where  the  one-way  wave  equation  method  fails, 
the  results  show  excellent  agreement.  Then  we  show  the  accuracy  of  the  method 
by  comparing  synthetic  seismograms  generated  by  this  method  with  those  gener¬ 
ated  by  a  finite  difference  algorithm  (Xie  and  Lay,  1994).  For  the  finite-difference 
method,  a  fourth-order  elastic  SH-wave  code  is  used  to  calculate  the  synthetic 
seismograms.  The  spatial  sampling  interval  is  0.25  km  and  the  time  interval  is 
0.025  second.  For  the  screen  method,  the  spatial  sampling  interval  is  0.25  km  in 
the  vertical  direction  and  the  screen  interval  is  1.0  km.  A  Gaussian  derivative  is 
used  as  the  source  time  function  for  both  methods.  Because  of  the  computational 
intensity  of  the  finite  difference  method,  we  did  the  comparison  at  short  propaga¬ 
tion  distances.  Shown  on  the  top  of  Figure  2  is  the  crustal  model  used  to  calculate 
synthetic  seismograms.  The  lower  part  shows  the  synthetic  seismograms  along  a 
vertical  profile  at  an  epicenter  distance  of  250  km.  The  thin  lines  are  from  the 
finite  difference  method  and  the  thick  lines  are  from  the  generalized  screen  method. 
The  source  is  located  at  a  depth  of  2  km.  Excellent  agreement  can  be  seen. 

Figure  3  shows  the  snap  shots  at  50  sec.  for  flat,  necking  and  broadening  crustal 
waveguides  (from  top  to  bottom,  respectively)  calculated  using  the  screen  method. 
The  source  is  located  at  the  top-left  corner  at  depth  2  km.  The  development  of  the 
mantle  wave  and  head  wave,  and  the  formation  of  crustal  guided  waves  as  multiple 
reflections  between  the  free  surface  and  Moho  discontinuity  can  be  clearly  seen. 
For  the  inhomogeneous  models,  wave  diffraction,  leakage  to  the  mantle,  wavefront 
distortion  and  an  increase  of  wavefield  complexity  can  be  seen  clearly  from  the 
snapshots.  It  can  be  seen  also  that  the  passage  through  a  narrow  crustal  segment 
has  greater  effect  on  Lg  leakage  than  the  passage  through  a  broad  segment.  In  the 
latter  case,  although  the  wavefronts  are  complicated  due  to  scattering  at  the  edges, 
most  of  the  energy  is  still  trapped  in  the  crust,  in  contrast  to  the  case  of  a  narrow 
passage,  in  which  a  large  percentage  of  energy  leaks  into  the  mantle.  This  example 


10 


demonstrates  the  potential  of  the  screen  method  as  a  tool  for  investigating  the 
path  effects  of  different  crustal  structures. 

The  following  example  shows  the  potential  capability  of  this  method  for  long 
distance  high-frequency  wave  propagation  in  a  laterally  varying  structure.  Figure 
4  shows  the  laterally  varying  crustal  model  used  in  the  calculation.  Figure  5  shows 
the  high-frequency  synthetic  seismograms  on  the  surface  at  distances  up  to  1000 
km  for  an  inhomogeneous  crustal  waveguide.  The  center  frequency  is  5  Hz  and 
the  maximum  frequency  is  10  Hz.  Shown  in  the  right  panel  are  synthetics  at  short 
distances  (up  to  500  km);  on  the  left  panel  are  synthetics  for  long  distances  (up 
to  1000  km).  In  this  case,  the  Lg  group  is  formed  by  multiple  reflections  by  the 
Moho  and  the  crustal  discontinuities.  This  example  demonstrates  the  potential 
capability  of  the  generalized  screen  methods.  In  comparison,  the  low-frequency  (/c 
=  1  Hz,  fmax  —  2  Hz)  synthetic  seismograms  are  shown  in  Figure  6.  It  is  clear 
that  without  high-frequency  content,  many  of  the  distinctive  features  associated 
with  the  Lg  measurements  cannot  be  adequately  modeled. 

2.5  The  Influence  of  Random  Heterogeneities  and  Rough 
Interfaces 

The  importance  of  small-scale  random  heterogeneities  to  seismic  wave  propagation 
is  well  known.  There  are  extensive  publications  on  this  subject  in  seismology.  How¬ 
ever,  the  role  of  random  heterogeneities  in  Lg  excitation,  propagation,  attenuation 
and  blockage,  is  still  unclear  due  to  the  complexity  of  the  problem.  The  theory 
of  wave  propagation  in  unbounded  random  media  has  been  well  developed.  How¬ 
ever,  for  waves  in  complex  crustal  waveguides  with  random  heterogeneities,  the 
theoretical  difficulties  are  overwhelming,  and  no  analytical  tools  are  available  for 
performing  realistic  calculations.  Numerical  simulation  is  an  attractive  alternative 
to  theory.  Some  finite-difference  simulations  have  been  conducted  (e.g.  Frankel 
and  Clayton,  1986;  Frankel,  1989;  Xie  and  Lay,  1994;  Jih,  1996).  Limited  by 
the  computation  power,  however,  the  finite-difference  results  are  often  for  short 
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Comparison  of  Uaue  Propagation 
in  Uarious  Crustal  Uaue  Guides 

C  t  =  50  sec) 


Figure  3:  Shown  from  the  top  to  the  bottom  are  snap-shots  at  50  sec.  for  flat, 
necking  and  broadening  crustal  waveguides.  The  development  of  mantle  wave  and 
head  wave,  and  the  formation  of  multiple  reflections  between  the  free  surface  and 
the  Moho  for  crustal  guided  waves  can  be  seen  clearly. 
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Parameters  of  Crustal  Model 


Layer 

Vs(km/sec) 

Density(g/cm ) 

Thickness(km) 

1 

3.00 

2.60 

1.00 

2 

3.46 

2.80 

14.00 

3 

3.76 

3.00 

22.00 

4 

4.65 

3.30 

Half-Space 

Crustal  Model 

Figure  4:  Flora-Asnes  crust  model.  Shown  in  the  upper  panel  are  model  parameters 
and  the  lower  panel  gives  the  geometry  of  the  model.  There  is  a  low  velocity- 
sedimentary  layer  in  the  top  1  km  and  a  velocity  discontinuity  at  a  depth  of  15 
km.  The  receivers  are  on  the  surface  and  shown  by  triangles. 
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Figure  5:  High-frequency  synthetic  seismograms  on  the  surface  at  distances  up  to 
1000  km  for  an  inhomogeneous  crustal  waveguide.  The  left  panel  shows  synthetics 
for  long  distances  (up  to  1000  km)  and  the  right  panel  shows  synthetics  of  short 
distances  (up  to  350  km). 
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Figure  6:  Low-frequency  synthetic  seismograms  on  the  surface  at  distances  up  to 
1000  km  for  an  inhomogeneous  crustal  waveguide  (Figure  4).  The  left  panel  shows 
synthetics  for  long  distances  (up  to  1000  km)  and  the  right  panel  shows  synthetics 
for  short  distances  (up  to  350  km). 


distances  or  low  frequencies.  Liu  and  Wu  (1994)  have  done  some  numerical  sim¬ 
ulation  using  the  phase-screen  method,  but  the  media  simulated  are  limited  to 
unbounded  media.  The  development  of  the  half-space  GSP  method  enables  us 
to  simulate  high-frequency  waves  propagating  in  complex  crustal  waveguides  to 
long  distances.  In  the  following,  we  will  show  two  examples  demonstrating  the 
capability  of  the  method. 

Figure  7  shows  a  heterogeneous  crustal  model  representing  a  ‘mountain  root’ 
with  small-scale  random  heterogeneities.  On  the  top  panel  is  the  velocity  model, 
and  the  comparisons  between  synthetic  seismograms  with  and  without  random 
heterogeneities  are  shown  on  the  middle  and  bottom  panels,  respectively.  The  het¬ 
erogeneities  have  an  exponential  correlation  function,  with  the  scale  length  ax  =  az 
=  1.6  km  (in  horizontal  and  vertical  directions,  respectively).  The  RMS  velocity 
perturbation  is  5%.  The  dominant  frequency  of  the  source  function  is  2  Hz.  Fig¬ 
ures  8  A  and  B  show  the  comparison  between  snapshots  for  waves  passing  through 
the  mountain  root’  with  and  without  random  heterogeneities,  respectively.  We 
see  that  random  heterogeneities  drastically  increase  the  leakage  of  waves  to  the 
mantle  and  the  complexity  of  the  waveforms.  Extensive  numerical  experiments 
will  be  conducted  to  study  the  different  influences  of  various  kinds  of  random  het¬ 
erogeneities.  It  has  been  shown  that  the  crustal  random  heterogeneities  are  highly 
anisotropic  in  scale  length  (Levander  and  Holliger,  1992;  Holliger  and  Levander, 
1992;  Wu  et  ah,  1994).  The  influences  of  the  random  heterogeneities  with  different 
stochastic  characteristics  can  be  explored  systematically. 

Figure  9  shows  the  comparison  of  synthetic  seismograms  and  snapshots  for 
models  with  and  without  a  rough  interface.  Shown  on  the  first  panel  is  a  crustal 
model  with  a  1  km  thick  low-velocity  top  layer.  The  bottom  of  the  low-velocity 
layer  is  a  rough  interface  with  0.2  km  RMS  random  depth  fluctuations.  The 
randomness  has  an  exponential  correlation  function  and  a  horizontal  scale  length 
of  0.5km.  The  source  is  located  at  2  km  depth  and  at  zero  distance,  the  receiver 
is  located  at  the  surface  and  at  a  distance  of  500  km.  The  second  and  third  panels 
show  the  comparison  of  synthetic  seismograms;  the  bottom  two  panels  show  the 
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Crust  mode 


comparison  of  snapshots.  We  see  that  rough  interfaces  of  sedimentary  layers  can 
also  increase  the  mantle  leakage  and  waveform  complexity  of  regional  waves. 

In  the  following,  we  will  test  and  discuss  how  these  small  scale  structures  affect 
Lg  energy  transport.  Before  applying  the  screen  propagator  method  to  a  random 
velocity  model,  we  first  check  the  validity  of  the  screen  method  by  comparing  its 
result  with  that  from  a  finite-difference  method.  The  upper  panel  of  Figure  10  gives 
a  random  velocity  model  with  5%  RMS  velocity  perturbation  and  a  correlation 
length  of  1  km  (exponential  correlation  function)  in  the  crust.  To  shorten  the 
computation  time,  we  used  a  16  km  thick  crust.  The  source  is  located  at  a  depth 
of  2  km  and  at  the  zero  distance.  A  receiver  array  is  put  on  the  free  surface  from 
zero  distance  to  250  km.  Both  finite-difference  method  and  screen  method  are 
used  to  calculate  the  synthetic  seismograms  for  the  same  model.  The  wave  energy 
are  then  calculated  from  these  synthetic  seismograms.  The  lower  panel  shows 
the  comparison  of  relative  energy  decay  curves  between  these  two  methods.  The 
solid  line  is  from  the  finite-difference  method  and  the  dotted  line  is  from  the  screen 
propagator  method.  The  agreement  between  these  two  methods  is  quite  well.  This 
proves  the  validity  of  the  screen  propagator  applied  to  the  energy  transfer  and 
partition  in  random  crustal  waveguides.  For  this  test  model,  the  FD  calculation 
has  Ax  =  Az  —  0.125fcm,  At  =  0.015  sec,  resulting  in  a  CPU  time  of  58  hours  on 
a  SUN  SPARC-4  work  station,  while  the  screen  method  has  Ax  =  Az  =  0.25 km, 
At  =  0.1  sec.,  and  a  CPU  time  of  0.5  hour  on  the  same  machine.  Both  calculations 
have  /0  =  1  Hz.  For  the  screen  method,  the  cutoff  frequency  fmax  =  2  Hz. 

The  next  test  is  for  a  normal  crustal  waveguide  with  a  5%  RMS  random  veloc¬ 
ity  perturbations  in  the  crust.  The  fluctuation  has  an  exponential  power  spectrum 
with  different  scales.  Figure  11  gives  the  attenuation  curves  for  different  charac¬ 
teristic  scales.  The  upper  panel  is  the  relative  total  energy,  which  is  the  energy 
contained  in  the  whole  seismogram  recorded  on  the  surface  versus  distance.  The 
solid  line  is  for  ka  =  1,  the  dotted  line  is  for  ka  =  10,  where  k  is  wavenumber 
and  a  is  the  correlation  length  of  the  random  perturbations;  and  the  dashed  line 
is  for  the  reference  model  without  heterogeneities.  We  see  that  for  the  reference 
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Figure  9:  Influence  of  a  rough  interface  on  Lg  propagation.  The  top  panel  is  a 
ciustal  model  with,  a  sedimentary  layer  for  which  the  bottom  is  a  rough  interface. 
The  source  and  receiver  are  located  at  (0,  2)  and  (500,  0).  The  second  and  third 
panels  show  the  comparison  of  synthetic  seismograms;  the  bottom  two  panels  show 
the  comparison  of  snapshots. 
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Crustal  model  with  exponential  random  media 
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Comparison  of  energy  attenuation  between  FD(solid  line)  and  GSP(dotted  line) 


Figure  10:  Comparison  of  energy  attenuation  curves  calculated  by  screen  propa¬ 
gator  method  and  finite-difference  method.  Shown  in  the  top  panel  is  the  velocity 
structure  including  random  heterogeneities  and  in  the  lower  panel  are  the  energy 
attenuation  curves. 
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model,  the  total  energy  remains  basically  constant  beyond  critical  distance,  which 
serves  as  a  checking  point  for  the  numerical  simulations.  The  lower  panel  gives  the 
logarithmic  relative  RMS  Lg  wave  amplitudes,  which  are  calculated  within  the  Lg 
window,  versus  distance.  Again,  the  solid,  dashed  and  dotted  lines  are  for  ka  =  1, 
ka  =  10  and  the  reference  model.  In  both  measurements  ka  =  1  lines  give  stronger 
attenuation  than  ka  =  10  lines.  We  see  also  that  the  Lg  amplitudes  attenuate 
more  rapidly  than  the  total  energy.  This  is  due  to  the  effect  of  scattering,  which 
diffuses  the  waves  out  of  the  Lg  window. 

The  existence  of  rough  surface  and  interfaces  has  effects  on  both  excitation 
and  propagation  of  Lg-waves.  First,  the  rough  surface  near  the  source  region  may 
affect  the  source  energy  partitioning.  Second,  along  the  entire  propagation  path, 
the  trapped  energy  in  the  waveguide  may  be  redirected  to  steeper  angles  due  to 
scattering,  causing  it  to  leak  into  the  upper  mantle.  This  can  cause  additional  Lg- 
wave  attenuation.  Similar  to  the  example  shown  above,  we  use  a  low  velocity  layer 
with  a  rough  lower  interface  to  simulate  the  effect  of  a  rough  surface.  In  Figure 
12,  the  upper  panel  is  the  velocity  structure  including  a  surface  sedimentary  layer 
with  a  random  interface.  The  sedimentary  layer  has  an  average  thickness  of  1 
km  and  a  RMS  depth  perturbation  of  0.2  km  with  correlation  length  of  1  km. 
The  middle  panel  is  the  energy  distribution  versus  distance  and  vertical  slowness. 

It  clearly  shows  that  with  the  existence  of  a  rough  interface,  considerable  energy 
moves  from  lower  vertical  slowness  to  higher  vertical  slowness.  In  other  words 
the  energy  propagation  directions  are  deflected  from  near  horizontal  to  steeper 
directions,  which  makes  more  energy  leak  into  the  upper  mantle  and  causes  extra 
Lg  wave  attenuation.  The  lower  panel  gives  the  energy  attenuation  curves  versus 
distance,  in  which  the  solid  line  is  for  the  model  with  a  rough  interface  and  the 

dotted  line  is  for  the  reference  model.  It  shows  clearly  the  effect  of  a  rough  interface 
on  Lg  attenuation. 
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Figure  11:  Attenuation  curves  for  a  flat  crust  with  random  heterogeneities  having 
different  characteristic  scales.  The  upper  panel  is  the  relative  energy  attenuation, 
and  the  lower  panel  is  the  logarithmic  relative  RMS  Lg  wave  amplitude  attenuation. 
The  solid  line  is  for  ka  =  1,  the  dotted  line  is  for  ka  =  10,  where  a  is  the  correlation 
length,  and  the  dashed  line  is  for  the  reference  model. 
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Angular  Spectrum  for  crustal  model  with  rough  sediment  interface 


lrV"ltelET  diSlribUti“  f°r  3  CIUStal  WaVeg“ide  m°del  sedimen- 

.  .  .  TOm  t<>P  *°  bottom>  crustal  structure;  energy  distribution  versus 
verf  ca  slowness  and  distance;  and  energy  attenuation  versus  d.stance.  Dotted 

me  m  .cates  propagat.on  through  reference  model;  solid  line  through  model  with 
rough  interface. 
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3  Global  Generalized  Reflection/Transmission 
Matrix  Method 

The  discretization  of  BIE  can  be  done  by  integration  of  the  Green’s  function  ei¬ 
ther  in  space  domain  (e.g.  Sanches-Sesma  and  Campillo,  1991),  or  in  wavenumber 
domain  using  the  discrete  wave  number  representation  (Bouchon,  1985;  Campillo 
and  Bouchon,  1985;  Chen,  1990,  1995,  1996).  In  the  latter  approach,  the  singu¬ 
larity  problem  of  the  Green’s  function  is  avoided  by  using  truncated  series.  The 
wavenumber  domain  BIE  has  another  advantage  that  it  can  be  easily  extended 
to  the  case  of  multilayered  media  with  irregular  interfaces.  In  Bouchon  et  al. 
(1989),  propagator  matrices  are  used  to  relate  equivalent  force  distributions  on 
neighboring  interfaces.  Chen  (1990,  1995,  1996)  related  the  fields  at  neighboring 
layers  by  global  reflection/transmission  coefficients  and  then  derived  the  global 
generalized  R/T  coefficients  to  relate  observations  and  sources.  In  these  methods, 
the  dimensionality  of  the  linear  system  to  solve  are  independent  of  the  number 
of  layers  involved.  The  computation  time  increases  only  linearly  with  the  num¬ 
ber  of  interfaces.  For  this  reason,  we  adopt  Chen’s  GGRTM  (Global  Generalized 
Reflection/Transmission  Matrix)  method  as  the  candidate  in  our  hybrid  method. 

The  GGRTM  can  be  viewed  as  an  extension  of  the  reflectivity  method  for  hori¬ 
zontally  layered  case  to  an  irregularly  layered  case,  and  it  has  been  demonstrated  to 
be  an  accurate  and  effective  method  to  simulate  seismic  waves  in  laterally  varying 
layered  media  (Chen,  1991,  1995,  1996).  For  example,  for  the  scattering  problem 
due  to  a  semi-circular  canyon  (shown  in  Figure  13),  GGRTM  can  provide  very  ac¬ 
curate  results.  Figures  14  and  15  show  the  comparisons  of  the  results  (solid  lines) 
computed  by  GGRTM  with  the  analytical  solutions  of  Trifunac  (1971,1973)  (dot¬ 
ted  lines)  for  various  normalized  frequencies,  showing  excellent  agreement  between 
them.  It  is  known  that  in  this  semi-circular  canyon  model,  there  are  two  sharp 
edges.  Many  other  methods,  e.g.,  Aki-Larner  method,  T-matrix  method  and  other 
high-frequency  asymptotic  methods,  fail  to  provide  correct  solutions. 
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Figure  13:  The  configuration  of  the  scattering  problem  due  to  a  semi-circular 
canyon  and  an  incident  plane  wave,  where  a  is  the  radius  of  the  canyon,  and  B  is 
the  angle  of  incident  wave. 

3.1  Connection  Formulation 

Assume  domain  II  is  the  model  space  we  are  interested  in  and  the  field  in  domain  I  is 
easy  to  calculate  by  other  less  expensive  methods.  According  to  the  representation 
theorem,  wave-fields  inside  domain  II  can  be  expressed  as 

+°o  ^ 

uH(x,u;)=  J  |r/(x')  +  u\x.')fj,(z')-^/'jGII(x.,x')dz'  (17) 

Where  u  and  t1  are  the  displacement  and  traction  fields  on  the  vertical  boundary 
surface  dividing  domain  I  and  II,  and  can  be  calculated  using  methods  valid  in 
domain  I,  \i  is  the  shear  rigidity,  and  Gn  is  the  Green’s  function  in  domain  II 
which  will  be  calculated  by  GGRTM. 


-2.  -2.  -1.  0.  t.  2.  3. 

distonce  (x/a) 


Figure  14:  The  frequency  responses  of  a  semi-circular  canyon  to  vertical  incident 
SH-wave  for  various  normalized  frequencies.  The  solid  lines  denote  our  results  and 
the  dotted  lines  denote  the  exact  solutions. 
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Figure  15:  The  same  response  as  Figure  14,  except  that  the  incident  angle  is  30°. 


3.2  Algorithm  of  Computing  Synthetic  Lg  Waves 

Having  the  connection  formulation,  we  can  use  GGRTM  to  compute  synthetic  Lg 
waves.  The  step-by-step  procedure  of  applying  GGRTM  to  computing  a  synthetic 
seismogram  in  a  general  irregularly  layered  medium  can  be  summarized  as  follows. 
Step  1 

Calculate  the  interface  matrices  for  each  interface  Q$  q{j)  Q(j)  qD  p(i) 
p(i)  p(i)  p(j)  .  r  .  ,  0  „  ,  .  ,  ’  iT’  U’  ^tt’  n’ 

W  ’  >  Ior  J- V,  •••,  W,  by  carrying  out  the  integrals  over  each  interface. 

These  interface  matrices  contain  the  structural  information  of  the  media  and  are 
defined  as  (Chen,  1990) 


_j  l/2 

(Q")"=2 Wl  /  +  (18) 

”  -LI  2 
_1  i/2 

(Q^)n  =  I  {e(i_1)(®)^-^)}exp[iSg(a:,n,m)]dx,  (19) 

—  L/2 
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_jU0'+1)I/(i+i) 

•  m 


2  fjtWv&L 


L/2 

j  { 1  +  [e0)(x)]2}12exp[fsW(x,  n,  m)}dx  , 


-L/2 


where  £W(z)  is  the  height  of  the  topography  for  the  jth  interface,  and 
K  =  2 nnL  U&  =  yfc'/3U))2  -  (fcn)2  ,  and  Im{i/W}  >  0, 

=  (km  ~  kn) x  +  u^\x)  -  +  r/i+1)  |a^)(x)|  , 


S*(X’ m)  =  ^  ~  k^x  +  -  &?)]  +  v®  |Afy-»)(*)| 


and 


S^)(af’ n’ m)  “ ^  “  **)*  -  x)(^)  -  eS;0] + ^  Ia^-^x)! 

for  j  =  1,  2,  ....  N. 

Where  Af^(x)  =  —  z^K 


Step  2 

®Iobal  modified  reflection  and/or  transmission  matrices,  {  R® 
Ttt  TU  Rti  }|  from  the  interface  matrices  using  the  following  formulas: 


R<?  T« 

T«  R» 


■Qii+,) 


p(i+i) 

'ru 


_q(-0  pO) 

q0+1)  pO'+l) 
'•Ht 


p{?  1  - 


EW 

max 


R«  =  -Qu  (Qff)"1  Ei'’„ . 

Where  Emin  and  E^x  are  diagonal  matrices  given  by 

E“  =  d, agonal  {exp[,uW(fW  -  {&•>)].  „  „  0>  ±1,  ±2, ...}  , 

and 


E-“  =  di“9onal  {exptiuW/^  _  ft,-D)];  „  =  ±1>  ±2>  j 

These  global  modified  reflection  and/or  transmission  matrices  describe  the  re¬ 
flection  and/or  transmission  effects  due  to  a  single  interface  regardless  of  the  in- 
fluences  from  other  existing  interfaces. 

Step  3 

a  ,  .)C°n^Ute  the  §Iobal  generalized  reflection  and/or  transmission  matrices,  T$ 

y  ’  u’  and  Rtt’  from  the  global  modified  reflection  and/or  transmission  ma¬ 
trices  through  following  recursive  formulas: 


and 


(28) 
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0 

[I  -  R<? RijHlr1T«')  ,  /or  j  =  N,  N  -  1 . 2,1.  (29) 

p(i)  i  'T'(i)-R(i+1)'-p(i) 

xvt4.  i-  J-U 

These  global  generalized  reflection  and/or  transmission  matrices  represent  the  total 
reflection  and/or  transmissions  due  to  the  multi-irregular  layers. 

Step  4 

Compute  the  expansion  coefficients  of  displacement  spectra  at  the  free  surface, 

2  =  £(0)(z), 


R 


(N+ 1)  _ 

It  ~ 


'T'O) 

1U 

pW 

■K-ti 


by  using  the  formula 
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fsd) 
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_i_  dd1)^2) 


+  T|1t)T(2) 
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(3) 


_L  i  rr(1)q^(2)  't(N)a(-'v+1)'1 

[  *  •  •  J  ^  A  •  •  •  X  I 


(30) 

where  is  the  equivalent  source  term  for  the  jth  layer  derived  by  the  represen¬ 
tation  theorem  and 


8»  =  {I  -  RfR!?}"  (sf  +  . 

L/2  ««>(*) 

(stj))n  =  ^~W7  S  dx  /  exp[— iknx  +  iu{J\z  -  &l)]dz  , 

“Vn  -L/2  {«-!)(*) 

L/2  &)(x) 

~  j  dx  J  fl3)(x,z)exp[-iknx-ii/^(z-^^)]dz  , 

-L/2  tti-V(x) 

for  j  =1,  2,  N,  N+l; 
and 

/(j)(x,  z)  =  |r(j)(x,  z)  +  x,  2)  j  • 

Step  5 


(31) 

(32) 


(33) 
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Calculate  the  displacement  spectra  at  the  free  surface  by  using  the  following 
formula: 

M 

W(o)[s,£(o)(x),a;]  =  E  «C)exp^m  +  zV(l1HAe(o)(a:)|}  .  (35) 

Taking  the  Fourier  transform  of  the  above  frequency  domain  solution,  we  can 
finally  obtain  the  time  domain  solution,  i.e.,  the  synthetic  seismogram. 

3.3  Numerical  Test 

To  test  the  validity  of  our  hybrid  method,  we  consider  a  trivial  case:  a  laterally 
homogeneous  layered  case.  This  problem  can  be  fully  solved  by  the  reflectivity 
method.  To  test  our  algorithm,  we  use  our  hybrid  method  to  synthesize  the  seis¬ 
mograms,  then  check  the  results  with  the  reflectivity  method.  The  test  model  is 
a  single  layer  crustal  model.  The  velocities  and  densities  of  the  crust  and  mantle 
are  3.5  km, /sec,  2.8  g/cm 3,  4.5  km /sec  and  3.2  g/cm3,  respectively.  The  thick¬ 
ness  of  the  crust  is  32  km ,  the  seismic  source  is  buried  at  zs=2  km  and  a:s=0 
km.  Receiver  is  placed  at  zo=0  km  and  xo=250  km.  The  connection  boundary  is 
located  at  ar=150  km.  The  synthetic  seismogram  of  the  reflectivity  method  is  plot¬ 
ted  in  Figure  16a.  The  synthetic  seismogram  of  GGRTM  is  shown  in  Figure  16b, 
Comparison  of  these  two  seismograms  shows  excellent  agreement  between  the  hy¬ 
brid  method  and  the  reflectivity  method,  confirming  the  validity  of  the  connection 
scheme  for  our  hybrid  method. 

The  computer  code  for  calculating  general  irregular  media  is  under  development 
at  this  stage,  and  is  expected  to  be  finished  soon.  We  will  then  calculate  synthetic 
Lg  waves  propagating  through  an  arbitrarily  irregular  layered  medium  to  study 
the  influence  of  surface  topography  and  interface  irregularities. 
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displacement  displacement 


Comparison  of  Synthetic  Seismograms 


Figure  16:  Comparison  of  synthetic  seismograms  for  a  laterally  homogeneous  lay¬ 
ered  crustal  model.  A:  synthetic  seismogram  from  reflectivity  method,  and  B: 
synthetic  seismogram  from  hybrid  method. 
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4 


A  Generalized  Screen  -  Boundary  Element 
Hybrid  Method 

4.1  Boundary  Element  Methods 

In  many  cases,  wave  propagation  involving  heterogeneous  media  can  be  formulated 
in  terms  of  boundary  integral  equations.  The  strict  boundary  methods  include  the 
direct  and  indirect  boundary  integral  equation  techniques.  The  direct  method  has 
been  widely  used  due  to  the  explicit  meaning  of  the  unknowns  in  the  formulation, 
while  the  indirect  method  formulates  problems  in  terms  of  force  or  force  moment 
boundary  densities.  These  two  methods  can  be'  implemented  either  in  the  space- 
time  domain  or  in  the  space-frequency  domain. 

For  many  applications,  an  alternative  approach  has  been  developed  by  the  com¬ 
bination  of  discrete  wavenumber  expansions  for  Green’s  functions  with  boundary 
integral  formulations  (Bouchon,  1985;  Campillo  and  Bouchon,  1985;  Bouchon  et 
al.,  1989).  In  this  case,  the  singularity  problem  of  the  Green’s  function  is  avoided 
by  using  wavenumber  integration.  This  technique  has  been  shown  to  be  accurate 
for  any  finite  frequency  and  can  be  extended  for  multilayered  media  with  irregular 
interfaces.  This  method  is  still  limited  to  the  case  of  low  frequencies  due  to  the 
computational  intensity. 

In  this  study,  we  select  the  direct  frequency-domain  BE  method  integrated 
with  the  screen  approach  to  model  the  effects  of  both  local  irregular  topography 
and  complex  crustal  structure  on  regional  wave  propagation.  In  application  to  the 
problem  of  regional  wave  propagation,  Gibson  and  Campillo  (1994)  have  used  the 
BE  method  for  frequencies  up  to  1  Hz.  The  computation  at  higher  frequencies  be¬ 
comes  extremely  time  consuming  because  matrices  with  large  size  must  be  inverted. 
In  fact,  the  propagation  properties  deduced  from  simulations  at  relatively  low  fre¬ 
quencies  on  the  order  of  1  Hz  show  very  different  characteristics  from  those  with 
higher  frequencies  (Wu  et  al.,  1996).  By  using  the  hybrid  scheme,  the  relatively 
short  sections  of  the  local  strong  heterogeneities  including  irregular  topography  in 
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the  crustal  waveguide  can  be  modeled  by  the  BE  method  to  high  frequencies,  and 
the  exterior  field  in  the  relative  weak  heterogeneous  media  of  large  volume  can  be 
calculated  by  the  screen  method. 

The  formulation  of  the  BE  method  can  be  briefly  described  as  follows.  Consider 
steady  state  scalar  wave  propagation  in  a  homogeneous  region  Cl  bounded  by  a 
boundary  T  .  Assuming  a  point  source  is  located  at  ro,  with  a  source  function 
S(uj)  ,  the  boundary  integral  equation  for  the  seismic  response  u( r)  at  location  r 
on  T  can  be  written  as 

C(r)tt(r)  +  J^u{r')— G(r,r’)dr  =  G(r,  r')— u(r')dr'  +  S(u)G(r,  r0),  (36) 

where  the  coefficients  C(r)  generally  depends  on  the  local  geometry  at  r  ,  G(r,  r') 
is  the  Green’s  function  for  the  homogeneous  region,  d/dn  denotes  differentiation 
with  respect  to  the  outward  normal  of  the  boundary  T,  and  S(u)  is  the  source 
spectrum. 

In  this  study  we  adopt  the  strict  frequency-space  domain  BE  method.  The 
boundary  T  is  discretized  into  a  finite  number  of  elements.  The  boundary  integral 
eq.  (1)  for  all  nodes  is  approximated  by  a  simultaneous  system  of  linear  equations. 
In  general,  the  coefficient  matrix  is  full-rank.  For  a  piecewise  homogeneous  media 
with  irregular  interfaces  (Fu  and  Mu,  1994;  Fu,  1996),  the  discretization  of  eq.  (1) 
can  be  done  in  each  subdomain,  and  then  all  equations  are  assembled  into  a  global 
matrix  equation  by  using  the  interface  conditions  of  continuity  for  displacements 
and  their  normal  gradients  across  all  interfaces.  This  global  matrix  is  sparse  or 
narrow-banded,  depending  on  the  structure  of  the  models.  Since  matrix  operations 
are  involved  and  the  matrix  for  each  frequency  component  must  be  inverted,  the 
BE  method  is  not  efficient  for  the  large  volume  problem.  This  problem  can  be 
circumvented  with  the  use  of  hybrid  methods. 

Figure  17  shows  the  application  of  the  BE  method  to  a  2-D  complex  salt  model 
in  Figure  17a.  The  medium  is  piecewise  homogeneous,  with  the  wave  velocities  in¬ 
dicated  in  the  figure.  The  dimensions  of  the  model  are  4000  m  horizontally  and  900 
m  vertically.  The  source  is  a  minimum-phase  wavelet  with  a  central  frequency  of 
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(b) 


Figure  17:  (a)  The  geometry  of  a  2-D  salt  model.  The  velocity  unit  is  m/s.  (b)  The  Synthetic 
acoustic  seismograms  calculated  with  the  BE  method. 
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20  Hz.  The  synthetic  seismograms  for  coincident  source-receiver  configurations  are 
displayed  in  Figure  17b  with  receivers  1  to  101  located  at  0  to  4000  m  along  the  x- 
axis.  Wave  fields  are  calculated  in  the  frequency  range  0  -  40  Hz  on  a  PC  computer 
with  Intel  166  MHz  Pentium  Processor.  A  variable  element  dimension  technique 
in  the  program  implementation  is  adopted  to  improve  the  computation  speed.  The 
element  dimension  for  each  frequency  is  computed  according  to  the  medium  ve¬ 
locity  and  the  frequency,  and  then  the  model  is  automatically  discretized.  This 
improves  the  efficiency  of  the  BE  method. 

4.2  A  GS-BE  Hybrid  Scheme 

Although  the  use  of  the  boundary  method  lags  far  behind  the  use  of  domain  meth¬ 
ods  (e.g.,  finite-difference  or  finite-element  methods)  in  engineering,  it  has  been 
extensively  used  to  study  the  effects  of  topography  or  sedimentary  basin  structures 
on  ground  motions  at  the  surface  during  the  last  two  decades  (Lee  and  Langston, 
1983;  Gaffet  and  Bouchon,  1989;  Sanchez-Sesma  et  ah,  1989;  Sanchez-Sesma  and 
Campillo,  1991;  1993;  Benites  and  Aki,  1989;  Papageorgiou  and  Kim,  1993;  Kim 
and  Papageorgiou,  1993;  and  has  gained  popularity  among  theoretical  seismolo¬ 
gists.  This  approach  has  also  been  used  to  study  the  Lg  blockage  problem  with 
limited  success.  Blockage  is  assumed  to  be  caused  by  coastlines,  mountains  and 
sudden  change  of  crustal  thickness.  However,  numerical  simulations  of  blockage 
by  large-scale  crustal  structures  by  using  finite  difference  and  boundary  methods 
have  not  succeeded  in  matching  the  observations  (Campillo.  1993;  Gibson  and 
Campillo,  1994).  Geological  structures  such  as  grabens  and  mountain  ranges  often 
lead  to  anomalous  attenuation,  even  extinction  of  the  Lg  phases.  However,  based 
on  the  synthetic  seismograms  using  these  methods,  the  Lg  waves  should  propagate 
with  a  much  smaller  degree  of  amplitude  reduction  than  is  observed.  One  of  the 
reasons  for  these  discrepancies  between  the  observed  and  modeling  results  is  that 
the  models  used  were  oversimplified  and  did  not  include  small-scale  irregularities. 

The  hybrid  method  will  combine  the  screen  method  and  the  BE  method  so 
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that  the  expensive  BE  calculations  can  be  limited  to  relatively  short  sections  with 
severe  surface  topography.  The  boundary  connection  technique  to  couple  the  fields 
calculated  by  the  screen  method  with  those  of  the  BE  method  will  play  a  crucial 
role  m  the  hybrid  method.  The  BE  method  is  a  well-tested  numerical  method. 
For  the  accuracy  of  the  screen  method,  extensive  numerical  tests  have  been  done 
(Liu  and  Wu,  1994;  Wu  et  ah,  1996;  Huang  and  Fehler,  1998)  in  comparison  with 
various  other  numerical  methods.  The  accuracy  of  the  hybrid  scheme  depends  on 
these  two  methods  and  their  connection. 


4.3  Boundary  Connection  Technique  for  the  Hybrid  Method 


The  problem  configuration  is  illustrated  in  Figure  18.  The  test  model  (Figure 
18a)  consists  of  an  irregular  free  surface  and  an  interface  as  a  single  layer  crustal 
model.  An  artificial  interface  F AB>  is  introduced  as  a  wavefield-connection  bound¬ 
ary  between  the  screen  and  BE  methods.  The  whole  model  is  divided  into  four 
subdomains,  flu  Fl2  in  the  crust,  and  f l[,  f V2  in  the  mantle.  The  left  boundaries 
of  flj  and  Fl[  ,  and  the  right  boundaries  of  ft2  and  Vt'2  are  assumed  to  extend  to 
infinity.  The  wavefields  in  and  Fl\  are  easy  to  calculate  by  the  screen  method. 
Its  output  field  ii0(r)  on  the  connection  interface  F  AB.  will  be  used  as  the  boundary 
condition  when  the  BE  method  is  used  to  calculate  wave  propagation  in  f}2  and  f 1' . 
The  output  field  is  received  along  the  next  connection  interface  FCD-,  and  will  be 
used  as  the  input  to  the  screen  method  in  the  next  propagation.  In  this  way,  a  long 
distance  propagation  of  the  Lg  phases  through  local  complex  waveguide  structures 
with  rough  topographic  features  can  be  simulated  by  the  hybrid  scheme. 

The  total  field  u(r)  on  the  boundaries  of  Q2  and  0'2  is  composed  of 


u(r)  = 


(37) 


u0(r)  +  u5(r) 

l  «,(r)  ,r£FAB.  ' 

The  boundary  integral  equation  for  fields  on  the  boundaries  of  fl2  can  be  obtained 
from 


C (r)u(r)  +  ^ u(r')— W 


(38) 
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The  surface  T  surrounding  subdomain  fl2  consists  of  the  connection  interface  Tab, 
the  upper  free  surface  IT,  the  lower  interface  T2,  and  the  right  boundary  Ycd- 
YcDi  in  computation,  is  assumed  to  be  transparent  and  can  be  handled  using  an 
infinite  element  absorbing  boundary  technique  (Fu  and  Wu,  1997).  Therefore,  the 
boundary  integrals  in  eq.  (38)  can  be  decomposed  into  integrals  over  segments  of 
the  boundary.  Considering  eq.  (37)  and  the  free  boundary  condition  on  Ti  ,  the 
integral  on  the  right  of  eq.  (38)  can  be  expressed  as 


In  order  to  solve  us(r)  and  dus(r)/dn  on  T2  ,  we  must  build  the  corresponding 
boundary  integral  equation  in  subdomain  ti'2  . 


C{r)us{r)  +  us(r')-^G(r,r')dr'  =  G{r,  r')-^us(r')dr' ,  (40) 

where  P  =  -f-  Y'2  +  YDD>  ,  and  a  sufficient  long  boundary  T BB>  is  used  with  its 

end  set  into  infinite  element  (Fu  and  Wu,  1997).  The  continuity  of  displacement 
and  its  normal  gradient  across  interface  Ybd  is  employed  when  eqs.  (38)  and 
(40)  are  combined  to  solve  the  problem.  To  solve  eq.  (38)  using  the  boundary 
condition,  the  normal  gradient  of  the  field  on  T^/  must  be  calculated.  It  is  easy 
to  calculate  du0(r)/dn  by  the  screen  method.  An  alternative  is  to  use  the  Rayleigh 
integral  representation  to  eliminate  the  term  du0(r)/dn.  In  addition,  since  there 
is  no  discontinuity  across  Y AB> ,  we  can  use  the  transparent  boundary  condition  on 
Y ab'  f°r  us(r).  By  solving  the  joint  boundary  integral  equations  of  02  and  Yl'2  ,  we 
can  obtain  the  wavefields  us(r)  on  Ti  and  T2  .  The  observation  field  along  YCD>  is 
calculated  explicitly  from  the  fields  on  the  boundaries. 

To  test  the  validity  of  the  connection  technique,  we  present  a  comparison  be¬ 
tween  the  wavefield  (Figure  18c)  obtained  using  the  BE  method  to  directly  cal¬ 
culate  wave  propagation  from  the  source  to  the  observation  surface  YCD>  and  the 
one  (Figure  18d)  by  the  hybrid  method  using  the  above  connection  scheme.  In 
both  cases,  the  source  wavelet  is  the  same,  with  the  dominant  frequency  at  20  Hz. 
First,  the  screen  method  is  used  to  compute  the  intermediate  wavefield  (Figure 
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Receiver  Number 


Figure  18:  Test  of  wavefield  connection,  (a)  Model  geometry  with  receivers  1-36  located  along 
TAB'  and  TCd’  respectively,  (b)  The  seismograms  computed  from  the  source  to  the  connection 
boundary  TAB,  .  (c)  The  seismograms  computed  from  the  source  directly  to  YCd<  ■  (d)  The 
seismograms  on  TCd>  computed  using  the  wavefields  on  as  the  incident  field.  Comparison 
of  wavefields  (c)  and  (d)  shows  the  validity  of  the  connection  technique. 
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18b)  on  V AB'  as  the  incident  field.  Then,  the  BE  method  is  used  to  calculate  wave 
propagation  from  T AB>  to  VCD'  .  The  dominant  arrivals  for  the  incident  field  at 
TAB>  consist  of  three  sets  of  waves.  The  first  arrival  is  the  direct  wave  propagating 
from  the  source  to  receivers.  The  second  is  the  reflected  wave  from  the  bottom 
interface,  and  the  third  is  the  reflected  wave  bounced  back  from  the  free  surface. 
From  the  seismograms  at  TCD >,  more  multiply  reflected  waves  between  the  free 
surface  and  the  interface  can  be  clearly  seen.  The  agreement  between  Figures  18c 
and  lSd  confirms  the  validity  of  the  connection  technique  for  the  hybrid  method. 

4.4  Numerical  Simulations  on  Surface  Topography 

In  this  section,  we  apply  the  hybrid  method  to  regional  wave  propagation  sim¬ 
ulations.  Figure  19  shows  a  laterally  varying  crust  model.  The  first  segment  is 
a  homogeneous  waveguide  200  km  in  length,  32  km  in  thickness,  with  a  SH  wave 
velocity  of  3.5  km/s  over  a  half-space  (4.5  km/s).  The  second  segment  is  a  complex 
waveguide  consisting  of  two  6-km  high,  20-km  long  mountain  ridges  centered  at 
240  km  and  290  km.  The  last  segment  represents  a  250  km  horizontal  waveguide 
of  32  km  thickness  with  small-scale  random  heterogeneities.  The  random  hetero¬ 
geneities  have  an  exponential  correlation  function  with  a  correlation  length  3  km 
in  both  horizontal  and  vertical  directions,  and  a  velocity  perturbation  6%.  The 
point  source  is  located  at  2  km  depth.  The  receivers  lie  along  both  the  free  surface 
and  vertical  profiles  at  different  distances  with  a  5  km  spacing.  Seismograms  are 
computed  in  the  frequency  range  0  -  2  Hz. 

First,  the  screen  method  is  used  to  compute  wave  propagation  from  the  source 
over  200  km  to  produce  an  initial  wavefield  on  the  first  connection  boundary  . 
The  synthetic  seismograms  are  shown  in  Figure  20  with  the  receivers  1-7  above  the 
Moho  and  8-10  under  the  Moho.  Multiple  reflections  within  the  crustal  waveguide 
can  be  seen  clearly  from  the  seismograms.  Then,  the  wavefields  on  are  used  as 
the  incident  fields,  and  the  BE  method  is  employed  to  calculate  wave  propagation 
for  the  next  120  km  through  the  segment  with  the  irregular  topography.  The 
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Figure  19.  A  laterally  varying  crustal  model.  The  shadowed  zone  represents  a  random  medium. 
The  receivers  are  indicated  by  small  triangles. 


variable  element  dimension  technique  with  6  elements  per  wavelength  is  used  in 
the  program  implementation  to  improve  the  computation  speed.  The  synthetic 
seismograms  shown  in  Figure  21  are  recorded  along  both  the  free  surface  and 
the  second  connection  boundary  TCd  at  320  km.  This  result  is  compared  to  the 
seismogram  in  Figure  22,  calculated  using  the  hybrid  method  for  the  same  segment 
of  the  model,  but  with  a  flat  free-surface.  Figure  23  shows  the  synthetic  seismogram 
obtained  using  only  the  screen  method  to  directly  calculate  wave  propagation  from 
the  source  to  the  connection  boundary  Tcd  for  the  model  with  the  flat  free-surface. 
The  agreement  between  Figures  22  and  17  shows  the  validity  of  the  hybrid  method 
for  a  long  distance  propagation. 

The  synthetic  seismograms  shown  in  Figure  21  are  more  complicated  than  those 
in  Figure  22.  The  effects  of  the  irregular  topography  on  the  guided  waves  can  be 
seen  from  the  comparison.  The  receivers  located  at  or  near  irregular  topographies 
have  anomalous  reflected  signals  due  to  the  focusing/defocusing  effects  and  the 
contamination  by  reverberations  within  these  topographic  structures  (local  site 
effects). 

The  most  characteristic  feature  of  the  effects  from  the  irregular  topography 
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Vertical  Seismograms  at  200  km 


Figure  20:  Synthetic  seismograms  along  the  first  connection  boundary  at  the  distance  of  200 
km  calculated  by  the  screen  method  to  produce  an  initial  wavefield. 

is  the  strong  scattering  by  the  topographic  structures.  This  scattering  has  two 
consequences  according  to  the  propagation  directions.  One  part  is  forescattering 
whose  active  time  window  is  close  to  that  of  the  wavetrain.  Therefore,  this  tends  to 
distort  the  wavetrain.  Another  part  is  backscattering  that  add  new  complexity  to 
the  waveguide  wavetrains  in  Figure  21.  In  the  surface  seismograms,  there  are  four 
groups  of  strong  arrivals  after  the  weak  head  waves:  they  are  the  direct,  primary, 
doubly  and  triply  reflected  waves  from  the  Moho.  We  see  very  little  backscattering 
for  the  direct  and  primary  reflected  waves,  but  see  clearly  the  backscattered  waves 
for  the  doubly  and  triply  reflected  waves.  Especially  for  the  triply  reflected  waves, 
the  backscatterings  around  the  two  mountain  ridges  at  240  km  and  290  km  are 
rather  strong  due  to  the  favorable  incident  angles  (steeper  angles)  in  this  case. 
From  the  vertical  profile  in  Figure  21  (bottom  panel)  we  see  more  mantle  waves 
compared  with  Figure  22  (without  topography).  This  is  due  to  the  leakage  of 
crustal  waves  into  the  mantle  by  the  mountain  ridge  scattering.  At  the  same 
time,  the  multiply  reflected  crustal  waves  have  been  weakened,  especially  for  the 


43 


Horizontal  Distance  (km) 


Surface  Seismograms 


Time  (s) 


Figure  22:  Same  as  Figure  21  except  the  free-surface  is  flat 


Vertical  Seismograms  at  570  km 


Figure  24:  Synthetic  seismoB„ms  along  a  vertical  profile  at  the  distance  of  570  km  calculated  by 
the  hybnd  method.  The  wavefields  Figure  21  is  used  as  the  incident  field  and  the  propagation 
in  the  random  medium  is  calculated  by  the  screen  method. 


quadruply  reflected  waves  (compared  to  Figure  22). 

In  summary,  the  scattering  by  local  irregular  topographies  leads  to  anomalous 
near-receiver  effects  and  tends  to  remove  some  energy  from  the  guided  waves, 
w  ich  causes  decay  of  amplitude  and  waveform  distortion.  These  scattered  waves 
partly  leak  to  the  mantle  and  partly  merge  into  the  crustal  guided  waves.  It  can 
be  expected  that  rough  surface  topography  with  scale  length  close  to  the  dominant 
wavelength  will  be  very  efficient  in  attenuating  Lg  waves. 

The  synthetic  seismograms  shown  in  Figure  24  are  obtained  by  using  the  wave- 
fields  on  VCD  in  Figure  21  as  the  input  to  the  screen  method  and  observing  along 
the  boundary  VEF  at  570  km.  The  importance  of  small-scale  random  hetero¬ 
geneities  to  seismic  wave  propagation  is  well  known.  The  role  of  random  hetero¬ 
geneities  in  Lg  wave  propagation  has  been  studied  by  Wu  et  al.(1996)  using  the 


half-space  generalized  screen  method,  which  shows  that  random  heterogeneities 
drastically  increase  the  leakage  of  waves  to  the  mantle  and  the  complexity  of  the 
waveforms. 

5  Conclusions 

Based  on  the  newly  developed  half-space  screen  propagator  method  formulated  by 
our  group,  hybrid  methods  have  been  proposed  and  tested  for  Lg  wave  simulation 
in  highly  complex  crustal  waveguides  including  surface  topography  for  the  two- 
dimensional  SH  case.  We  have  derived  the  connection  formulas  for  our  hybrid 
schemes  and  their  validity  has  been  proved  by  numerical  tests.  Generalized  screen 
method  and  both  boundary  integral  equation  and  boundary  element  methods  have 
been  tested  for  the  waveguide  environment.  The  excellent  agreement  between 
seismograms  for  direct  propagation  and  propagation  using  the  connection  formulas 
proves  the  correctness  of  the  theory  and  the  connection  formulations.  Numerical 
simulations  on  the  influence  of  surface  topography,  sedimentary  layers  with  rough 
bottom,  and  small-scale  random  heterogeneities  demonstrate  the  feasibility  of  the 
methodology.  It  is  also  shown  that  rough  surface  topography  and  the  rough  bottom 
of  sedimentary  layers  with  scale  close  to  the  dominant  wavelength  can  efficiently 
attenuate  Lg  waves.  The  next  step  is  to  develop  the  corresponding  theory  and 
algorithms  for  2D  P-SV  and  3D  elastic  wave  problems. 
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